home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / hyperg_1F1.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  58.9 KB  |  2,007 lines
  1. /* specfunc/hyperg_1F1.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_elementary.h>
  26. #include <gsl/gsl_sf_exp.h>
  27. #include <gsl/gsl_sf_bessel.h>
  28. #include <gsl/gsl_sf_gamma.h>
  29. #include <gsl/gsl_sf_laguerre.h>
  30. #include <gsl/gsl_sf_hyperg.h>
  31.  
  32. #include "error.h"
  33. #include "hyperg.h"
  34.  
  35. #define _1F1_INT_THRESHOLD (100.0*GSL_DBL_EPSILON)
  36.  
  37.  
  38. /* Asymptotic result for 1F1(a, b, x)  x -> -Infinity.
  39.  * Assumes b-a != neg integer and b != neg integer.
  40.  */
  41. static
  42. int
  43. hyperg_1F1_asymp_negx(const double a, const double b, const double x,
  44.                      gsl_sf_result * result)
  45. {
  46.   gsl_sf_result lg_b;
  47.   gsl_sf_result lg_bma;
  48.   double sgn_b;
  49.   double sgn_bma;
  50.  
  51.   int stat_b   = gsl_sf_lngamma_sgn_e(b,   &lg_b,   &sgn_b);
  52.   int stat_bma = gsl_sf_lngamma_sgn_e(b-a, &lg_bma, &sgn_bma);
  53.  
  54.   if(stat_b == GSL_SUCCESS && stat_bma == GSL_SUCCESS) {
  55.     gsl_sf_result F;
  56.     int stat_F = gsl_sf_hyperg_2F0_series_e(a, 1.0+a-b, -1.0/x, -1, &F);
  57.     if(F.val != 0) {
  58.       double ln_term_val = a*log(-x);
  59.       double ln_term_err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + fabs(ln_term_val));
  60.       double ln_pre_val = lg_b.val - lg_bma.val - ln_term_val;
  61.       double ln_pre_err = lg_b.err + lg_bma.err + ln_term_err;
  62.       int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
  63.                                             sgn_bma*sgn_b*F.val, F.err,
  64.                                             result);
  65.       return GSL_ERROR_SELECT_2(stat_e, stat_F);
  66.     }
  67.     else {
  68.       result->val = 0.0;
  69.       result->err = 0.0;
  70.       return stat_F;
  71.     }
  72.   }
  73.   else {
  74.     DOMAIN_ERROR(result);
  75.   }
  76. }
  77.  
  78.  
  79. /* Asymptotic result for 1F1(a, b, x)  x -> +Infinity
  80.  * Assumes b != neg integer and a != neg integer
  81.  */
  82. static
  83. int
  84. hyperg_1F1_asymp_posx(const double a, const double b, const double x,
  85.                       gsl_sf_result * result)
  86. {
  87.   gsl_sf_result lg_b;
  88.   gsl_sf_result lg_a;
  89.   double sgn_b;
  90.   double sgn_a;
  91.  
  92.   int stat_b = gsl_sf_lngamma_sgn_e(b, &lg_b, &sgn_b);
  93.   int stat_a = gsl_sf_lngamma_sgn_e(a, &lg_a, &sgn_a);
  94.  
  95.   if(stat_a == GSL_SUCCESS && stat_b == GSL_SUCCESS) {
  96.     gsl_sf_result F;
  97.     int stat_F = gsl_sf_hyperg_2F0_series_e(b-a, 1.0-a, 1.0/x, -1, &F);
  98.     if(stat_F == GSL_SUCCESS && F.val != 0) {
  99.       double lnx = log(x);
  100.       double ln_term_val = (a-b)*lnx;
  101.       double ln_term_err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + fabs(b)) * fabs(lnx)
  102.                          + 2.0 * GSL_DBL_EPSILON * fabs(a-b);
  103.       double ln_pre_val = lg_b.val - lg_a.val + ln_term_val + x;
  104.       double ln_pre_err = lg_b.err + lg_a.err + ln_term_err + 2.0 * GSL_DBL_EPSILON * fabs(x);
  105.       int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
  106.                                             sgn_a*sgn_b*F.val, F.err,
  107.                                             result);
  108.       return GSL_ERROR_SELECT_2(stat_e, stat_F);
  109.     }
  110.     else {
  111.       result->val = 0.0;
  112.       result->err = 0.0;
  113.       return stat_F;
  114.     }
  115.   }
  116.   else {
  117.     DOMAIN_ERROR(result);
  118.   }
  119. }
  120.  
  121.  
  122. /* Asymptotic result for x < 2b-4a, 2b-4a large.
  123.  * [Abramowitz+Stegun, 13.5.21]
  124.  *
  125.  * assumes 0 <= x/(2b-4a) <= 1
  126.  */
  127. static
  128. int
  129. hyperg_1F1_large2bm4a(const double a, const double b, const double x, gsl_sf_result * result)
  130. {
  131.   double eta    = 2.0*b - 4.0*a;
  132.   double cos2th = x/eta;
  133.   double sin2th = 1.0 - cos2th;
  134.   double th = acos(sqrt(cos2th));
  135.   double pre_h  = 0.25*M_PI*M_PI*eta*eta*cos2th*sin2th;
  136.   gsl_sf_result lg_b;
  137.   int stat_lg = gsl_sf_lngamma_e(b, &lg_b);
  138.   double t1 = 0.5*(1.0-b)*log(0.25*x*eta);
  139.   double t2 = 0.25*log(pre_h);
  140.   double lnpre_val = lg_b.val + 0.5*x + t1 - t2;
  141.   double lnpre_err = lg_b.err + 2.0 * GSL_DBL_EPSILON * (fabs(0.5*x) + fabs(t1) + fabs(t2));
  142.   double s1 = sin(a*M_PI);
  143.   double s2 = sin(0.25*eta*(2.0*th - sin(2.0*th)) + 0.25*M_PI);
  144.   double ser_val = s1 + s2;
  145.   double ser_err = 2.0 * GSL_DBL_EPSILON * (fabs(s1) + fabs(s2));
  146.   int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
  147.                                         ser_val, ser_err,
  148.                                         result);
  149.   return GSL_ERROR_SELECT_2(stat_e, stat_lg);
  150. }
  151.  
  152.  
  153. /* Luke's rational approximation.
  154.  * See [Luke, Algorithms for the Computation of Mathematical Functions, p.182]
  155.  *
  156.  * Like the case of the 2F1 rational approximations, these are
  157.  * probably guaranteed to converge for x < 0, barring gross
  158.  * numerical instability in the pre-asymptotic regime.
  159.  */
  160. static
  161. int
  162. hyperg_1F1_luke(const double a, const double c, const double xin,
  163.                 gsl_sf_result * result)
  164. {
  165.   const double RECUR_BIG = 1.0e+50;
  166.   const int nmax = 5000;
  167.   int n = 3;
  168.   const double x  = -xin;
  169.   const double x3 = x*x*x;
  170.   const double t0 = a/c;
  171.   const double t1 = (a+1.0)/(2.0*c);
  172.   const double t2 = (a+2.0)/(2.0*(c+1.0));
  173.   double F = 1.0;
  174.   double prec;
  175.  
  176.   double Bnm3 = 1.0;                                  /* B0 */
  177.   double Bnm2 = 1.0 + t1 * x;                         /* B1 */
  178.   double Bnm1 = 1.0 + t2 * x * (1.0 + t1/3.0 * x);    /* B2 */
  179.  
  180.   double Anm3 = 1.0;                                                      /* A0 */
  181.   double Anm2 = Bnm2 - t0 * x;                                            /* A1 */
  182.   double Anm1 = Bnm1 - t0*(1.0 + t2*x)*x + t0 * t1 * (c/(c+1.0)) * x*x;   /* A2 */
  183.  
  184.   while(1) {
  185.     double npam1 = n + a - 1;
  186.     double npcm1 = n + c - 1;
  187.     double npam2 = n + a - 2;
  188.     double npcm2 = n + c - 2;
  189.     double tnm1  = 2*n - 1;
  190.     double tnm3  = 2*n - 3;
  191.     double tnm5  = 2*n - 5;
  192.     double F1 =  (n-a-2) / (2*tnm3*npcm1);
  193.     double F2 =  (n+a)*npam1 / (4*tnm1*tnm3*npcm2*npcm1);
  194.     double F3 = -npam2*npam1*(n-a-2) / (8*tnm3*tnm3*tnm5*(n+c-3)*npcm2*npcm1);
  195.     double E  = -npam1*(n-c-1) / (2*tnm3*npcm2*npcm1);
  196.  
  197.     double An = (1.0+F1*x)*Anm1 + (E + F2*x)*x*Anm2 + F3*x3*Anm3;
  198.     double Bn = (1.0+F1*x)*Bnm1 + (E + F2*x)*x*Bnm2 + F3*x3*Bnm3;
  199.     double r = An/Bn;
  200.  
  201.     prec = fabs((F - r)/F);
  202.     F = r;
  203.  
  204.     if(prec < GSL_DBL_EPSILON || n > nmax) break;
  205.  
  206.     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
  207.       An   /= RECUR_BIG;
  208.       Bn   /= RECUR_BIG;
  209.       Anm1 /= RECUR_BIG;
  210.       Bnm1 /= RECUR_BIG;
  211.       Anm2 /= RECUR_BIG;
  212.       Bnm2 /= RECUR_BIG;
  213.       Anm3 /= RECUR_BIG;
  214.       Bnm3 /= RECUR_BIG;
  215.     }
  216.     else if(fabs(An) < 1.0/RECUR_BIG || fabs(Bn) < 1.0/RECUR_BIG) {
  217.       An   *= RECUR_BIG;
  218.       Bn   *= RECUR_BIG;
  219.       Anm1 *= RECUR_BIG;
  220.       Bnm1 *= RECUR_BIG;
  221.       Anm2 *= RECUR_BIG;
  222.       Bnm2 *= RECUR_BIG;
  223.       Anm3 *= RECUR_BIG;
  224.       Bnm3 *= RECUR_BIG;
  225.     }
  226.  
  227.     n++;
  228.     Bnm3 = Bnm2;
  229.     Bnm2 = Bnm1;
  230.     Bnm1 = Bn;
  231.     Anm3 = Anm2;
  232.     Anm2 = Anm1;
  233.     Anm1 = An;
  234.   }
  235.  
  236.   result->val  = F;
  237.   result->err  = 2.0 * fabs(F * prec);
  238.   result->err += 2.0 * GSL_DBL_EPSILON * (n-1.0) * fabs(F);
  239.  
  240.   return GSL_SUCCESS;
  241. }
  242.  
  243.  
  244. /* Series for 1F1(1,b,x)
  245.  * b > 0
  246.  */
  247. static
  248. int
  249. hyperg_1F1_1_series(const double b, const double x, gsl_sf_result * result)
  250. {
  251.   double sum_val = 1.0;
  252.   double sum_err = 0.0;
  253.   double term = 1.0;
  254.   double n    = 1.0;
  255.   while(fabs(term/sum_val) > 2.0*GSL_DBL_EPSILON) {
  256.     term *= x/(b+n-1);
  257.     sum_val += term;
  258.     sum_err += 2.0 * 4.0 * GSL_DBL_EPSILON * fabs(term);
  259.     n += 1.0;
  260.   }
  261.   result->val  = sum_val;
  262.   result->err  = sum_err;
  263.   result->err += 2.0 * fabs(term);
  264.   return GSL_SUCCESS;
  265. }
  266.  
  267.  
  268. /* 1F1(1,b,x)
  269.  * b >= 1, b integer
  270.  */
  271. static
  272. int
  273. hyperg_1F1_1_int(const int b, const double x, gsl_sf_result * result)
  274. {
  275.   if(b < 1) {
  276.     DOMAIN_ERROR(result);
  277.   }
  278.   else if(b == 1) {
  279.     return gsl_sf_exp_e(x, result);
  280.   }
  281.   else if(b == 2) {
  282.     return gsl_sf_exprel_e(x, result);
  283.   }
  284.   else if(b == 3) {
  285.     return gsl_sf_exprel_2_e(x, result);
  286.   }
  287.   else {
  288.     return gsl_sf_exprel_n_e(b-1, x, result);
  289.   }
  290. }
  291.  
  292.  
  293. /* 1F1(1,b,x)
  294.  * b >=1, b real
  295.  *
  296.  * checked OK: [GJ] Thu Oct  1 16:46:35 MDT 1998
  297.  */
  298. static
  299. int
  300. hyperg_1F1_1(const double b, const double x, gsl_sf_result * result)
  301. {
  302.   double ax = fabs(x);
  303.   double ib = floor(b + 0.1);
  304.  
  305.   if(b < 1.0) {
  306.     DOMAIN_ERROR(result);
  307.   }
  308.   else if(b == 1.0) {
  309.     return gsl_sf_exp_e(x, result);
  310.   }
  311.   else if(b >= 1.4*ax) {
  312.     return hyperg_1F1_1_series(b, x, result);
  313.   }
  314.   else if(fabs(b - ib) < _1F1_INT_THRESHOLD && ib < INT_MAX) {
  315.     return hyperg_1F1_1_int((int)ib, x, result);
  316.   }
  317.   else if(x > 0.0) {
  318.     if(x > 100.0 && b < 0.75*x) {
  319.       return hyperg_1F1_asymp_posx(1.0, b, x, result);
  320.     }
  321.     else if(b < 1.0e+05) {
  322.       /* Recurse backward on b, from a
  323.        * chosen offset point. For x > 0,
  324.        * which holds here, this should
  325.        * be a stable direction.
  326.        */
  327.       const double off = ceil(1.4*x-b) + 1.0;
  328.       double bp = b + off;
  329.       gsl_sf_result M;
  330.       int stat_s = hyperg_1F1_1_series(bp, x, &M);
  331.       const double err_rat = M.err / fabs(M.val);
  332.       while(bp > b+0.1) {
  333.         /* M(1,b-1) = x/(b-1) M(1,b) + 1 */
  334.         bp -= 1.0;
  335.         M.val  = 1.0 + x/bp * M.val;
  336.       }
  337.       result->val  = M.val;
  338.       result->err  = err_rat * fabs(M.val);
  339.       result->err += 2.0 * GSL_DBL_EPSILON * (fabs(off)+1.0) * fabs(M.val);
  340.       return stat_s;
  341.     }
  342.     else {
  343.       return hyperg_1F1_large2bm4a(1.0, b, x, result);
  344.     }
  345.   }
  346.   else {
  347.     /* x <= 0 and b not large compared to |x|
  348.      */
  349.     if(ax < 10.0 && b < 10.0) {
  350.       return hyperg_1F1_1_series(b, x, result);
  351.     }
  352.     else if(ax >= 100.0 && GSL_MAX_DBL(fabs(2.0-b),1.0) < 0.99*ax) {
  353.       return hyperg_1F1_asymp_negx(1.0, b, x, result);
  354.     }
  355.     else {
  356.       return hyperg_1F1_luke(1.0, b, x, result);
  357.     }
  358.   }
  359. }
  360.  
  361.  
  362. /* 1F1(a,b,x)/Gamma(b) for b->0
  363.  * [limit of Abramowitz+Stegun 13.3.7]
  364.  */
  365. static
  366. int
  367. hyperg_1F1_renorm_b0(const double a, const double x, gsl_sf_result * result)
  368. {
  369.   double eta = a*x;
  370.   if(eta > 0.0) {
  371.     double root_eta = sqrt(eta);
  372.     gsl_sf_result I1_scaled;
  373.     int stat_I = gsl_sf_bessel_I1_scaled_e(2.0*root_eta, &I1_scaled);
  374.     if(I1_scaled.val <= 0.0) {
  375.       result->val = 0.0;
  376.       result->err = 0.0;
  377.       return GSL_ERROR_SELECT_2(stat_I, GSL_EDOM);
  378.     }
  379.     else {
  380.       const double lnr_val = 0.5*x + 0.5*log(eta) + fabs(x) + log(I1_scaled.val);
  381.       const double lnr_err = GSL_DBL_EPSILON * (1.5*fabs(x) + 1.0) + fabs(I1_scaled.err/I1_scaled.val);
  382.       return gsl_sf_exp_err_e(lnr_val, lnr_err, result);
  383.     }
  384.   }
  385.   else if(eta == 0.0) {
  386.     result->val = 0.0;
  387.     result->err = 0.0;
  388.     return GSL_SUCCESS;
  389.   }
  390.   else {
  391.     /* eta < 0 */
  392.     double root_eta = sqrt(-eta);
  393.     gsl_sf_result J1;
  394.     int stat_J = gsl_sf_bessel_J1_e(2.0*root_eta, &J1);
  395.     if(J1.val <= 0.0) {
  396.       result->val = 0.0;
  397.       result->err = 0.0;
  398.       return GSL_ERROR_SELECT_2(stat_J, GSL_EDOM);
  399.     }
  400.     else {
  401.       const double t1 = 0.5*x;
  402.       const double t2 = 0.5*log(-eta);
  403.       const double t3 = fabs(x);
  404.       const double t4 = log(J1.val);
  405.       const double lnr_val = t1 + t2 + t3 + t4;
  406.       const double lnr_err = GSL_DBL_EPSILON * (1.5*fabs(x) + 1.0) + fabs(J1.err/J1.val);
  407.       gsl_sf_result ex;
  408.       int stat_e = gsl_sf_exp_err_e(lnr_val, lnr_err, &ex);
  409.       result->val = -ex.val;
  410.       result->err =  ex.err;
  411.       return stat_e;
  412.     }
  413.   }
  414.   
  415. }
  416.  
  417.  
  418. /* 1F1'(a,b,x)/1F1(a,b,x)
  419.  * Uses Gautschi's version of the CF.
  420.  * [Gautschi, Math. Comp. 31, 994 (1977)]
  421.  *
  422.  * Supposedly this suffers from the "anomalous convergence"
  423.  * problem when b < x. I have seen anomalous convergence
  424.  * in several of the continued fractions associated with
  425.  * 1F1(a,b,x). This particular CF formulation seems stable
  426.  * for b > x. However, it does display a painful artifact
  427.  * of the anomalous convergence; the convergence plateaus
  428.  * unless b >>> x. For example, even for b=1000, x=1, this
  429.  * method locks onto a ratio which is only good to about
  430.  * 4 digits. Apparently the rest of the digits are hiding
  431.  * way out on the plateau, but finite-precision lossage
  432.  * means you will never get them.
  433.  */
  434. #if 0
  435. static
  436. int
  437. hyperg_1F1_CF1_p(const double a, const double b, const double x, double * result)
  438. {
  439.   const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  440.   const int maxiter = 5000;
  441.   int n = 1;
  442.   double Anm2 = 1.0;
  443.   double Bnm2 = 0.0;
  444.   double Anm1 = 0.0;
  445.   double Bnm1 = 1.0;
  446.   double a1 = 1.0;
  447.   double b1 = 1.0;
  448.   double An = b1*Anm1 + a1*Anm2;
  449.   double Bn = b1*Bnm1 + a1*Bnm2;
  450.   double an, bn;
  451.   double fn = An/Bn;
  452.  
  453.   while(n < maxiter) {
  454.     double old_fn;
  455.     double del;
  456.     n++;
  457.     Anm2 = Anm1;
  458.     Bnm2 = Bnm1;
  459.     Anm1 = An;
  460.     Bnm1 = Bn;
  461.     an = (a+n)*x/((b-x+n-1)*(b-x+n));
  462.     bn = 1.0;
  463.     An = bn*Anm1 + an*Anm2;
  464.     Bn = bn*Bnm1 + an*Bnm2;
  465.  
  466.     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
  467.       An /= RECUR_BIG;
  468.       Bn /= RECUR_BIG;
  469.       Anm1 /= RECUR_BIG;
  470.       Bnm1 /= RECUR_BIG;
  471.       Anm2 /= RECUR_BIG;
  472.       Bnm2 /= RECUR_BIG;
  473.     }
  474.  
  475.     old_fn = fn;
  476.     fn = An/Bn;
  477.     del = old_fn/fn;
  478.     
  479.     if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
  480.   }
  481.  
  482.   *result = a/(b-x) * fn;
  483.  
  484.   if(n == maxiter)
  485.     GSL_ERROR ("error", GSL_EMAXITER);
  486.   else
  487.     return GSL_SUCCESS;
  488. }
  489. #endif /* 0 */
  490.  
  491.  
  492. /* 1F1'(a,b,x)/1F1(a,b,x)
  493.  * Uses Gautschi's series transformation of the
  494.  * continued fraction. This is apparently the best
  495.  * method for getting this ratio in the stable region.
  496.  * The convergence is monotone and supergeometric
  497.  * when b > x.
  498.  * Assumes a >= -1.
  499.  */
  500. static
  501. int
  502. hyperg_1F1_CF1_p_ser(const double a, const double b, const double x, double * result)
  503. {
  504.   if(a == 0.0) {
  505.     *result = 0.0;
  506.     return GSL_SUCCESS;
  507.   }
  508.   else {
  509.     const int maxiter = 5000;
  510.     double sum  = 1.0;
  511.     double pk   = 1.0;
  512.     double rhok = 0.0;
  513.     int k;
  514.     for(k=1; k<maxiter; k++) {
  515.       double ak = (a + k)*x/((b-x+k-1.0)*(b-x+k));
  516.       rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0+rhok));
  517.       pk  *= rhok;
  518.       sum += pk;
  519.       if(fabs(pk/sum) < 2.0*GSL_DBL_EPSILON) break;
  520.     }
  521.     *result = a/(b-x) * sum;
  522.     if(k == maxiter)
  523.       GSL_ERROR ("error", GSL_EMAXITER);
  524.     else
  525.       return GSL_SUCCESS;
  526.   }
  527. }
  528.  
  529.  
  530. /* 1F1(a+1,b,x)/1F1(a,b,x)
  531.  *
  532.  * I think this suffers from typical "anomalous convergence".
  533.  * I could not find a region where it was truly useful.
  534.  */
  535. #if 0
  536. static
  537. int
  538. hyperg_1F1_CF1(const double a, const double b, const double x, double * result)
  539. {
  540.   const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  541.   const int maxiter = 5000;
  542.   int n = 1;
  543.   double Anm2 = 1.0;
  544.   double Bnm2 = 0.0;
  545.   double Anm1 = 0.0;
  546.   double Bnm1 = 1.0;
  547.   double a1 = b - a - 1.0;
  548.   double b1 = b - x - 2.0*(a+1.0);
  549.   double An = b1*Anm1 + a1*Anm2;
  550.   double Bn = b1*Bnm1 + a1*Bnm2;
  551.   double an, bn;
  552.   double fn = An/Bn;
  553.  
  554.   while(n < maxiter) {
  555.     double old_fn;
  556.     double del;
  557.     n++;
  558.     Anm2 = Anm1;
  559.     Bnm2 = Bnm1;
  560.     Anm1 = An;
  561.     Bnm1 = Bn;
  562.     an = (a + n - 1.0) * (b - a - n);
  563.     bn = b - x - 2.0*(a+n);
  564.     An = bn*Anm1 + an*Anm2;
  565.     Bn = bn*Bnm1 + an*Bnm2;
  566.  
  567.     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
  568.       An /= RECUR_BIG;
  569.       Bn /= RECUR_BIG;
  570.       Anm1 /= RECUR_BIG;
  571.       Bnm1 /= RECUR_BIG;
  572.       Anm2 /= RECUR_BIG;
  573.       Bnm2 /= RECUR_BIG;
  574.     }
  575.  
  576.     old_fn = fn;
  577.     fn = An/Bn;
  578.     del = old_fn/fn;
  579.     
  580.     if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
  581.   }
  582.  
  583.   *result = fn;
  584.   if(n == maxiter)
  585.     GSL_ERROR ("error", GSL_EMAXITER);
  586.   else
  587.     return GSL_SUCCESS;
  588. }
  589. #endif /* 0 */
  590.  
  591.  
  592. /* 1F1(a,b+1,x)/1F1(a,b,x)
  593.  *
  594.  * This seemed to suffer from "anomalous convergence".
  595.  * However, I have no theory for this recurrence.
  596.  */
  597. #if 0
  598. static
  599. int
  600. hyperg_1F1_CF1_b(const double a, const double b, const double x, double * result)
  601. {
  602.   const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  603.   const int maxiter = 5000;
  604.   int n = 1;
  605.   double Anm2 = 1.0;
  606.   double Bnm2 = 0.0;
  607.   double Anm1 = 0.0;
  608.   double Bnm1 = 1.0;
  609.   double a1 = b + 1.0;
  610.   double b1 = (b + 1.0) * (b - x);
  611.   double An = b1*Anm1 + a1*Anm2;
  612.   double Bn = b1*Bnm1 + a1*Bnm2;
  613.   double an, bn;
  614.   double fn = An/Bn;
  615.  
  616.   while(n < maxiter) {
  617.     double old_fn;
  618.     double del;
  619.     n++;
  620.     Anm2 = Anm1;
  621.     Bnm2 = Bnm1;
  622.     Anm1 = An;
  623.     Bnm1 = Bn;
  624.     an = (b + n) * (b + n - 1.0 - a) * x;
  625.     bn = (b + n) * (b + n - 1.0 - x);
  626.     An = bn*Anm1 + an*Anm2;
  627.     Bn = bn*Bnm1 + an*Bnm2;
  628.  
  629.     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
  630.       An /= RECUR_BIG;
  631.       Bn /= RECUR_BIG;
  632.       Anm1 /= RECUR_BIG;
  633.       Bnm1 /= RECUR_BIG;
  634.       Anm2 /= RECUR_BIG;
  635.       Bnm2 /= RECUR_BIG;
  636.     }
  637.  
  638.     old_fn = fn;
  639.     fn = An/Bn;
  640.     del = old_fn/fn;
  641.     
  642.     if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
  643.   }
  644.  
  645.   *result = fn;
  646.   if(n == maxiter)
  647.     GSL_ERROR ("error", GSL_EMAXITER);
  648.   else
  649.     return GSL_SUCCESS;
  650. }
  651. #endif /* 0 */
  652.  
  653.  
  654. /* 1F1(a,b,x)
  655.  * |a| <= 1, b > 0
  656.  */
  657. static
  658. int
  659. hyperg_1F1_small_a_bgt0(const double a, const double b, const double x, gsl_sf_result * result)
  660. {
  661.   const double bma = b-a;
  662.   const double oma = 1.0-a;
  663.   const double ap1mb = 1.0+a-b;
  664.   const double abs_bma = fabs(bma);
  665.   const double abs_oma = fabs(oma);
  666.   const double abs_ap1mb = fabs(ap1mb);
  667.  
  668.   const double ax = fabs(x);
  669.  
  670.   if(a == 0.0) {
  671.     result->val = 1.0;
  672.     result->err = 0.0;
  673.     return GSL_SUCCESS;
  674.   }
  675.   else if(a == 1.0 && b >= 1.0) {
  676.     return hyperg_1F1_1(b, x, result);
  677.   }
  678.   else if(a == -1.0) {
  679.     result->val  = 1.0 + a/b * x;
  680.     result->err  = GSL_DBL_EPSILON * (1.0 + fabs(a/b * x));
  681.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  682.     return GSL_SUCCESS;
  683.   }
  684.   else if(b >= 1.4*ax) {
  685.     return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  686.   }
  687.   else if(x > 0.0) {
  688.     if(x > 100.0 && abs_bma*abs_oma < 0.5*x) {
  689.       return hyperg_1F1_asymp_posx(a, b, x, result);
  690.     }
  691.     else if(b < 5.0e+06) {
  692.       /* Recurse backward on b from
  693.        * a suitably high point.
  694.        */
  695.       const double b_del = ceil(1.4*x-b) + 1.0;
  696.       double bp = b + b_del;
  697.       gsl_sf_result r_Mbp1;
  698.       gsl_sf_result r_Mb;
  699.       double Mbp1;
  700.       double Mb;
  701.       double Mbm1;
  702.       int stat_0 = gsl_sf_hyperg_1F1_series_e(a, bp+1.0, x, &r_Mbp1);
  703.       int stat_1 = gsl_sf_hyperg_1F1_series_e(a, bp,     x, &r_Mb);
  704.       const double err_rat = fabs(r_Mbp1.err/r_Mbp1.val) + fabs(r_Mb.err/r_Mb.val);
  705.       Mbp1 = r_Mbp1.val;
  706.       Mb   = r_Mb.val;
  707.       while(bp > b+0.1) {
  708.         /* Do backward recursion. */
  709.         Mbm1 = ((x+bp-1.0)*Mb - x*(bp-a)/bp*Mbp1)/(bp-1.0);
  710.         bp -= 1.0;
  711.     Mbp1 = Mb;
  712.     Mb   = Mbm1;
  713.       }
  714.       result->val  = Mb;
  715.       result->err  = err_rat * (fabs(b_del)+1.0) * fabs(Mb);
  716.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mb);
  717.       return GSL_ERROR_SELECT_2(stat_0, stat_1);
  718.     }
  719.     else {
  720.       return hyperg_1F1_large2bm4a(a, b, x, result);
  721.     }
  722.   }
  723.   else {
  724.     /* x < 0 and b not large compared to |x|
  725.      */
  726.     if(ax < 10.0 && b < 10.0) {
  727.       return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  728.     }
  729.     else if(ax >= 100.0 && GSL_MAX(abs_ap1mb,1.0) < 0.99*ax) {
  730.       return hyperg_1F1_asymp_negx(a, b, x, result);
  731.     }
  732.     else {
  733.       return hyperg_1F1_luke(a, b, x, result);
  734.     }
  735.   }
  736. }
  737.  
  738.  
  739. /* 1F1(b+eps,b,x)
  740.  * |eps|<=1, b > 0
  741.  */
  742. static
  743. int
  744. hyperg_1F1_beps_bgt0(const double eps, const double b, const double x, gsl_sf_result * result)
  745. {
  746.   if(b > fabs(x) && fabs(eps) < GSL_SQRT_DBL_EPSILON) {
  747.     /* If b-a is very small and x/b is not too large we can
  748.      * use this explicit approximation.
  749.      *
  750.      * 1F1(b+eps,b,x) = exp(ax/b) (1 - eps x^2 (v2 + v3 x + ...) + ...)
  751.      *
  752.      *   v2 = a/(2b^2(b+1))
  753.      *   v3 = a(b-2a)/(3b^3(b+1)(b+2))
  754.      *   ...
  755.      *
  756.      * See [Luke, Mathematical Functions and Their Approximations, p.292]
  757.      *
  758.      * This cannot be used for b near a negative integer or zero.
  759.      * Also, if x/b is large the deviation from exp(x) behaviour grows.
  760.      */
  761.     double a = b + eps;
  762.     gsl_sf_result exab;
  763.     int stat_e = gsl_sf_exp_e(a*x/b, &exab);
  764.     double v2 = a/(2.0*b*b*(b+1.0));
  765.     double v3 = a*(b-2.0*a)/(3.0*b*b*b*(b+1.0)*(b+2.0));
  766.     double v  = v2 + v3 * x;
  767.     double f  = (1.0 - eps*x*x*v);
  768.     result->val  = exab.val * f;
  769.     result->err  = exab.err * fabs(f);
  770.     result->err += fabs(exab.val) * GSL_DBL_EPSILON * (1.0 + fabs(eps*x*x*v));
  771.     result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val);
  772.     return stat_e;
  773.   }
  774.   else {
  775.     /* Otherwise use a Kummer transformation to reduce
  776.      * it to the small a case.
  777.      */
  778.     gsl_sf_result Kummer_1F1;
  779.     int stat_K = hyperg_1F1_small_a_bgt0(-eps, b, -x, &Kummer_1F1);
  780.     if(Kummer_1F1.val != 0.0) {
  781.       int stat_e = gsl_sf_exp_mult_err_e(x, 2.0*GSL_DBL_EPSILON*fabs(x),
  782.                                             Kummer_1F1.val, Kummer_1F1.err,
  783.                                             result);
  784.       return GSL_ERROR_SELECT_2(stat_e, stat_K);
  785.     }
  786.     else {
  787.       result->val = 0.0;
  788.       result->err = 0.0;
  789.       return stat_K;
  790.     }
  791.   }
  792. }
  793.  
  794.  
  795. /* 1F1(a,2a,x) = Gamma(a + 1/2) E(x) (|x|/4)^(-a+1/2) scaled_I(a-1/2,|x|/2)
  796.  *
  797.  * E(x) = exp(x) x > 0
  798.  *      = 1      x < 0
  799.  *
  800.  * a >= 1/2
  801.  */
  802. static
  803. int
  804. hyperg_1F1_beq2a_pos(const double a, const double x, gsl_sf_result * result)
  805. {
  806.   if(x == 0.0) {
  807.     result->val = 1.0;
  808.     result->err = 0.0;
  809.     return GSL_SUCCESS;
  810.   }
  811.   else {
  812.     gsl_sf_result I;
  813.     int stat_I = gsl_sf_bessel_Inu_scaled_e(a-0.5, 0.5*fabs(x), &I);
  814.     gsl_sf_result lg;
  815.     int stat_g = gsl_sf_lngamma_e(a + 0.5, &lg);
  816.     double ln_term   = (0.5-a)*log(0.25*fabs(x));
  817.     double lnpre_val = lg.val + GSL_MAX_DBL(x,0.0) + ln_term;
  818.     double lnpre_err = lg.err + GSL_DBL_EPSILON * (fabs(ln_term) + fabs(x));
  819.     int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
  820.                           I.val, I.err,
  821.                           result);
  822.     return GSL_ERROR_SELECT_3(stat_e, stat_g, stat_I);
  823.   }
  824. }
  825.  
  826.  
  827. /* Determine middle parts of diagonal recursion along b=2a
  828.  * from two endpoints, i.e.
  829.  *
  830.  * given:  M(a,b)      and  M(a+1,b+2)
  831.  * get:    M(a+1,b+1)  and  M(a,b+1)
  832.  */
  833. #if 0
  834. inline
  835. static
  836. int
  837. hyperg_1F1_diag_step(const double a, const double b, const double x,
  838.                      const double Mab, const double Map1bp2,
  839.                      double * Map1bp1, double * Mabp1)
  840. {
  841.   if(a == b) {
  842.     *Map1bp1 = Mab;
  843.     *Mabp1   = Mab - x/(b+1.0) * Map1bp2;
  844.   }
  845.   else {
  846.     *Map1bp1 = Mab - x * (a-b)/(b*(b+1.0)) * Map1bp2;
  847.     *Mabp1   = (a * *Map1bp1 - b * Mab)/(a-b);
  848.   }
  849.   return GSL_SUCCESS;
  850. }
  851. #endif /* 0 */
  852.  
  853.  
  854. /* Determine endpoint of diagonal recursion.
  855.  *
  856.  * given:  M(a,b)    and  M(a+1,b+2)
  857.  * get:    M(a+1,b)  and  M(a+1,b+1)
  858.  */
  859. #if 0
  860. inline
  861. static
  862. int
  863. hyperg_1F1_diag_end_step(const double a, const double b, const double x,
  864.                          const double Mab, const double Map1bp2,
  865.                          double * Map1b, double * Map1bp1)
  866. {
  867.   *Map1bp1 = Mab - x * (a-b)/(b*(b+1.0)) * Map1bp2;
  868.   *Map1b   = Mab + x/b * *Map1bp1;
  869.   return GSL_SUCCESS;
  870. }
  871. #endif /* 0 */
  872.  
  873.  
  874. /* Handle the case of a and b both positive integers.
  875.  * Assumes a > 0 and b > 0.
  876.  */
  877. static
  878. int
  879. hyperg_1F1_ab_posint(const int a, const int b, const double x, gsl_sf_result * result)
  880. {
  881.   double ax = fabs(x);
  882.  
  883.   if(a == b) {
  884.     return gsl_sf_exp_e(x, result);             /* 1F1(a,a,x) */
  885.   }
  886.   else if(a == 1) {
  887.     return gsl_sf_exprel_n_e(b-1, x, result);   /* 1F1(1,b,x) */
  888.   }
  889.   else if(b == a + 1) {
  890.     gsl_sf_result K;
  891.     int stat_K = gsl_sf_exprel_n_e(a, -x, &K);  /* 1F1(1,1+a,-x) */
  892.     int stat_e = gsl_sf_exp_mult_err_e(x, 2.0 * GSL_DBL_EPSILON * fabs(x),
  893.                           K.val, K.err,
  894.                           result);
  895.     return GSL_ERROR_SELECT_2(stat_e, stat_K);
  896.   }
  897.   else if(a == b + 1) {
  898.     gsl_sf_result ex;
  899.     int stat_e = gsl_sf_exp_e(x, &ex);
  900.     result->val  = ex.val * (1.0 + x/b);
  901.     result->err  = ex.err * (1.0 + x/b);
  902.     result->err += ex.val * GSL_DBL_EPSILON * (1.0 + fabs(x/b));
  903.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  904.     return stat_e;
  905.   }
  906.   else if(a == b + 2) {
  907.     gsl_sf_result ex;
  908.     int stat_e = gsl_sf_exp_e(x, &ex);
  909.     double poly  = (1.0 + x/b*(2.0 + x/(b+1.0)));
  910.     result->val  = ex.val * poly;
  911.     result->err  = ex.err * fabs(poly);
  912.     result->err += ex.val * GSL_DBL_EPSILON * (1.0 + fabs(x/b) * (2.0 + fabs(x/(b+1.0))));
  913.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  914.     return stat_e;
  915.   }
  916.   else if(b == 2*a) {
  917.     return hyperg_1F1_beq2a_pos(a, x, result);  /* 1F1(a,2a,x) */
  918.   }
  919.   else if(   ( b < 10 && a < 10 && ax < 5.0 )
  920.           || ( b > a*ax )
  921.       || ( b > a && ax < 5.0 )
  922.     ) {
  923.     return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  924.   }
  925.   else if(b > a && b >= 2*a + x) {
  926.     /* Use the Gautschi CF series, then
  927.      * recurse backward to a=0 for normalization.
  928.      * This will work for either sign of x.
  929.      */
  930.     double rap;
  931.     int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
  932.     double ra = 1.0 + x/a * rap;
  933.     double Ma    = GSL_SQRT_DBL_MIN;
  934.     double Map1 = ra * Ma;
  935.     double Mnp1 = Map1;
  936.     double Mn    = Ma;
  937.     double Mnm1;
  938.     int n;
  939.     for(n=a; n>0; n--) {
  940.       Mnm1 = (n * Mnp1 - (2*n-b+x) * Mn) / (b-n);
  941.       Mnp1 = Mn;
  942.       Mn   = Mnm1;
  943.     }
  944.     result->val = Ma/Mn;
  945.     result->err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + 1.0) * fabs(Ma/Mn);
  946.     return stat_CF1;
  947.   }
  948.   else if(b > a && b < 2*a + x && b > x) {
  949.     /* Use the Gautschi series representation of
  950.      * the continued fraction. Then recurse forward
  951.      * to the a=b line for normalization. This will
  952.      * work for either sign of x, although we do need
  953.      * to check for b > x, for when x is positive.
  954.      */
  955.     double rap;
  956.     int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
  957.     double ra = 1.0 + x/a * rap;
  958.     gsl_sf_result ex;
  959.     int stat_ex;
  960.  
  961.     double Ma   = GSL_SQRT_DBL_MIN;
  962.     double Map1 = ra * Ma;
  963.     double Mnm1 = Ma;
  964.     double Mn    = Map1;
  965.     double Mnp1;
  966.     int n;
  967.     for(n=a+1; n<b; n++) {
  968.       Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
  969.       Mnm1 = Mn;
  970.       Mn   = Mnp1;
  971.     }
  972.  
  973.     stat_ex = gsl_sf_exp_e(x, &ex);  /* 1F1(b,b,x) */
  974.     result->val  = ex.val * Ma/Mn;
  975.     result->err  = ex.err * fabs(Ma/Mn);
  976.     result->err += 4.0 * GSL_DBL_EPSILON * (fabs(b-a)+1.0) * fabs(result->val);
  977.     return GSL_ERROR_SELECT_2(stat_ex, stat_CF1);
  978.   }
  979.   else if(x >= 0.0) {
  980.  
  981.     if(b < a) {
  982.       /* The point b,b is below the b=2a+x line.
  983.        * Forward recursion on a from b,b+1 is possible.
  984.        * Note that a > b + 1 as well, since we already tried a = b + 1.
  985.        */
  986.       if(x + log(fabs(x/b)) < GSL_LOG_DBL_MAX-2.0) {
  987.         double ex = exp(x);
  988.         int n;
  989.         double Mnm1 = ex;          /* 1F1(b,b,x)   */
  990.         double Mn   = ex * (1.0 + x/b);   /* 1F1(b+1,b,x) */
  991.         double Mnp1;
  992.         for(n=b+1; n<a; n++) {
  993.           Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
  994.           Mnm1 = Mn;
  995.           Mn   = Mnp1;
  996.         }
  997.         result->val  = Mn;
  998.     result->err  = (x + 1.0) * GSL_DBL_EPSILON * fabs(Mn);
  999.     result->err *= fabs(a-b)+1.0;
  1000.         return GSL_SUCCESS;
  1001.       }
  1002.       else {
  1003.         OVERFLOW_ERROR(result);
  1004.       }
  1005.     }
  1006.     else {
  1007.       /* b > a
  1008.        * b < 2a + x 
  1009.        * b <= x (otherwise we would have finished above)
  1010.        *
  1011.        * Gautschi anomalous convergence region. However, we can
  1012.        * recurse forward all the way from a=0,1 because we are
  1013.        * always underneath the b=2a+x line.
  1014.        */
  1015.       gsl_sf_result r_Mn;
  1016.       double Mnm1 = 1.0;    /* 1F1(0,b,x) */
  1017.       double Mn;        /* 1F1(1,b,x)  */
  1018.       double Mnp1;
  1019.       int n;
  1020.       gsl_sf_exprel_n_e(b-1, x, &r_Mn);
  1021.       Mn = r_Mn.val;
  1022.       for(n=1; n<a; n++) {
  1023.         Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
  1024.          Mnm1 = Mn;
  1025.          Mn   = Mnp1;
  1026.       }
  1027.       result->val  = Mn;
  1028.       result->err  = fabs(Mn) * (1.0 + fabs(a)) * fabs(r_Mn.err / r_Mn.val);
  1029.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn);
  1030.       return GSL_SUCCESS;
  1031.     }
  1032.   }
  1033.   else {
  1034.     /* x < 0
  1035.      * b < a (otherwise we would have tripped one of the above)
  1036.      */
  1037.  
  1038.     if(a <= 0.5*(b-x) || a >= -x) {
  1039.       /* Gautschi continued fraction is in the anomalous region,
  1040.        * so we must find another way. We recurse down in b,
  1041.        * from the a=b line.
  1042.        */
  1043.       double ex = exp(x);
  1044.       double Manp1 = ex;
  1045.       double Man   = ex * (1.0 + x/(a-1.0));
  1046.       double Manm1;
  1047.       int n;
  1048.       for(n=a-1; n>b; n--) {
  1049.         Manm1 = (-n*(1-n-x)*Man - x*(n-a)*Manp1)/(n*(n-1.0));
  1050.         Manp1 = Man;
  1051.         Man = Manm1;
  1052.       }
  1053.       result->val  = Man;
  1054.       result->err  = (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(Man);
  1055.       result->err *= fabs(b-a)+1.0;
  1056.       return GSL_SUCCESS;
  1057.     }
  1058.     else {
  1059.       /* Pick a0 such that b ~= 2a0 + x, then
  1060.        * recurse down in b from a0,a0 to determine
  1061.        * the values near the line b=2a+x. Then recurse
  1062.        * forward on a from a0.
  1063.        */
  1064.       int a0 = ceil(0.5*(b-x));
  1065.       double Ma0b;    /* M(a0,b)   */
  1066.       double Ma0bp1;  /* M(a0,b+1) */
  1067.       double Ma0p1b;  /* M(a0+1,b) */
  1068.       double Mnm1;
  1069.       double Mn;
  1070.       double Mnp1;
  1071.       int n;
  1072.       {
  1073.         double ex = exp(x);
  1074.         double Ma0np1 = ex;
  1075.         double Ma0n   = ex * (1.0 + x/(a0-1.0));
  1076.         double Ma0nm1;
  1077.         for(n=a0-1; n>b; n--) {
  1078.           Ma0nm1 = (-n*(1-n-x)*Ma0n - x*(n-a0)*Ma0np1)/(n*(n-1.0));
  1079.           Ma0np1 = Ma0n;
  1080.           Ma0n = Ma0nm1;
  1081.         }
  1082.     Ma0bp1 = Ma0np1;
  1083.         Ma0b   = Ma0n;
  1084.     Ma0p1b = (b*(a0+x)*Ma0b + x*(a0-b)*Ma0bp1)/(a0*b);
  1085.       }
  1086.  
  1087.       /* Initialise the recurrence correctly BJG */
  1088.  
  1089.       if (a0 >= a)
  1090.         { 
  1091.           Mn = Ma0b;
  1092.         }
  1093.       else if (a0 + 1>= a)
  1094.         {
  1095.           Mn = Ma0p1b;
  1096.         }
  1097.       else
  1098.         {
  1099.           Mnm1 = Ma0b;
  1100.           Mn   = Ma0p1b;
  1101.  
  1102.           for(n=a0+1; n<a; n++) {
  1103.             Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
  1104.             Mnm1 = Mn;
  1105.             Mn   = Mnp1;
  1106.           }
  1107.         }
  1108.  
  1109.       result->val  = Mn;
  1110.       result->err  = (fabs(x) + 1.0) * GSL_DBL_EPSILON *  fabs(Mn);
  1111.       result->err *= fabs(b-a)+1.0;
  1112.       return GSL_SUCCESS;
  1113.     }
  1114.   }
  1115. }
  1116.  
  1117.  
  1118. /* Evaluate a <= 0, a integer, cases directly. (Polynomial; Horner)
  1119.  * When the terms are all positive, this
  1120.  * must work. We will assume this here.
  1121.  */
  1122. static
  1123. int
  1124. hyperg_1F1_a_negint_poly(const int a, const double b, const double x, gsl_sf_result * result)
  1125. {
  1126.   if(a == 0) {
  1127.     result->val = 1.0;
  1128.     result->err = 1.0;
  1129.     return GSL_SUCCESS;
  1130.   }
  1131.   else {
  1132.     int N = -a;
  1133.     double poly = 1.0;
  1134.     int k;
  1135.     for(k=N-1; k>=0; k--) {
  1136.       double t = (a+k)/(b+k) * (x/(k+1));
  1137.       double r = t + 1.0/poly;
  1138.       if(r > 0.9*GSL_DBL_MAX/poly) {
  1139.         OVERFLOW_ERROR(result);
  1140.       }
  1141.       else {
  1142.         poly *= r;  /* P_n = 1 + t_n P_{n-1} */
  1143.       }
  1144.     }
  1145.     result->val = poly;
  1146.     result->err = 2.0 * (sqrt(N) + 1.0) * GSL_DBL_EPSILON * fabs(poly);
  1147.     return GSL_SUCCESS;
  1148.   }
  1149. }
  1150.  
  1151.  
  1152. /* Evaluate negative integer a case by relation
  1153.  * to Laguerre polynomials. This is more general than
  1154.  * the direct polynomial evaluation, but is safe
  1155.  * for all values of x.
  1156.  *
  1157.  * 1F1(-n,b,x) = n!/(b)_n Laguerre[n,b-1,x]
  1158.  *           = n B(b,n) Laguerre[n,b-1,x]
  1159.  *
  1160.  * assumes b is not a negative integer
  1161.  */
  1162. static
  1163. int
  1164. hyperg_1F1_a_negint_lag(const int a, const double b, const double x, gsl_sf_result * result)
  1165. {
  1166.   const int n = -a;
  1167.  
  1168.   gsl_sf_result lag;
  1169.   const int stat_l = gsl_sf_laguerre_n_e(n, b-1.0, x, &lag);
  1170.   if(b < 0.0) {
  1171.     gsl_sf_result lnfact;
  1172.     gsl_sf_result lng1;
  1173.     gsl_sf_result lng2;
  1174.     double s1, s2;
  1175.     const int stat_f  = gsl_sf_lnfact_e(n, &lnfact);
  1176.     const int stat_g1 = gsl_sf_lngamma_sgn_e(b + n, &lng1, &s1);
  1177.     const int stat_g2 = gsl_sf_lngamma_sgn_e(b, &lng2, &s2);
  1178.     const double lnpre_val = lnfact.val - (lng1.val - lng2.val);
  1179.     const double lnpre_err = lnfact.err + lng1.err + lng2.err
  1180.       + 2.0 * GSL_DBL_EPSILON * fabs(lnpre_val);
  1181.     const int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
  1182.                                                 s1*s2*lag.val, lag.err,
  1183.                                                 result);
  1184.     return GSL_ERROR_SELECT_5(stat_e, stat_l, stat_g1, stat_g2, stat_f);
  1185.   }
  1186.   else {
  1187.     gsl_sf_result lnbeta;
  1188.     gsl_sf_lnbeta_e(b, n, &lnbeta);
  1189.     if(fabs(lnbeta.val) < 0.1) {
  1190.       /* As we have noted, when B(x,y) is near 1,
  1191.        * evaluating log(B(x,y)) is not accurate.
  1192.        * Instead we evaluate B(x,y) directly.
  1193.        */
  1194.       const double ln_term_val = log(1.25*n);
  1195.       const double ln_term_err = 2.0 * GSL_DBL_EPSILON * ln_term_val;
  1196.       gsl_sf_result beta;
  1197.       int stat_b = gsl_sf_beta_e(b, n, &beta);
  1198.       int stat_e = gsl_sf_exp_mult_err_e(ln_term_val, ln_term_err,
  1199.                                             lag.val, lag.err,
  1200.                                             result);
  1201.       result->val *= beta.val/1.25;
  1202.       result->err *= beta.val/1.25;
  1203.       return GSL_ERROR_SELECT_3(stat_e, stat_l, stat_b);
  1204.     }
  1205.     else {
  1206.       /* B(x,y) was not near 1, so it is safe to use
  1207.        * the logarithmic values.
  1208.        */
  1209.       const double ln_n = log(n);
  1210.       const double ln_term_val = lnbeta.val + ln_n;
  1211.       const double ln_term_err = lnbeta.err + 2.0 * GSL_DBL_EPSILON * fabs(ln_n);
  1212.       int stat_e = gsl_sf_exp_mult_err_e(ln_term_val, ln_term_err,
  1213.                                             lag.val, lag.err,
  1214.                                             result);
  1215.       return GSL_ERROR_SELECT_2(stat_e, stat_l);
  1216.     }
  1217.   }
  1218. }
  1219.  
  1220.  
  1221. /* Handle negative integer a case for x > 0 and
  1222.  * generic b.
  1223.  *
  1224.  * Combine [Abramowitz+Stegun, 13.6.9 + 13.6.27]
  1225.  * M(-n,b,x) = (-1)^n / (b)_n U(-n,b,x) = n! / (b)_n Laguerre^(b-1)_n(x)
  1226.  */
  1227. #if 0
  1228. static
  1229. int
  1230. hyperg_1F1_a_negint_U(const int a, const double b, const double x, gsl_sf_result * result)
  1231. {
  1232.   const int n = -a;
  1233.   const double sgn = ( GSL_IS_ODD(n) ? -1.0 : 1.0 );
  1234.   double sgpoch;
  1235.   gsl_sf_result lnpoch;
  1236.   gsl_sf_result U;
  1237.   const int stat_p = gsl_sf_lnpoch_sgn_e(b, n, &lnpoch, &sgpoch);
  1238.   const int stat_U = gsl_sf_hyperg_U_e(-n, b, x, &U);
  1239.   const int stat_e = gsl_sf_exp_mult_err_e(-lnpoch.val, lnpoch.err,
  1240.                             sgn * sgpoch * U.val, U.err,
  1241.                             result);
  1242.   return GSL_ERROR_SELECT_3(stat_e, stat_U, stat_p);
  1243. }
  1244. #endif
  1245.  
  1246.  
  1247. /* Assumes a <= -1,  b <= -1, and b <= a.
  1248.  */
  1249. static
  1250. int
  1251. hyperg_1F1_ab_negint(const int a, const int b, const double x, gsl_sf_result * result)
  1252. {
  1253.   if(x == 0.0) {
  1254.     result->val = 1.0;
  1255.     result->err = 0.0;
  1256.     return GSL_SUCCESS;
  1257.   }
  1258.   else if(x > 0.0) {
  1259.     return hyperg_1F1_a_negint_poly(a, b, x, result);
  1260.   }
  1261.   else {
  1262.     /* Apply a Kummer transformation to make x > 0 so
  1263.      * we can evaluate the polynomial safely. Of course,
  1264.      * this assumes b <= a, which must be true for
  1265.      * a<0 and b<0, since otherwise the thing is undefined.
  1266.      */
  1267.     gsl_sf_result K;
  1268.     int stat_K = hyperg_1F1_a_negint_poly(b-a, b, -x, &K);
  1269.     int stat_e = gsl_sf_exp_mult_err_e(x, 2.0 * GSL_DBL_EPSILON * fabs(x),
  1270.                           K.val, K.err,
  1271.                           result);
  1272.     return GSL_ERROR_SELECT_2(stat_e, stat_K);
  1273.   }
  1274. }
  1275.  
  1276.  
  1277. /* [Abramowitz+Stegun, 13.1.3]
  1278.  *
  1279.  * M(a,b,x) = Gamma(1+a-b)/Gamma(2-b) x^(1-b) *
  1280.  *            { Gamma(b)/Gamma(a) M(1+a-b,2-b,x) - (b-1) U(1+a-b,2-b,x) }
  1281.  *
  1282.  * b not an integer >= 2
  1283.  * a-b not a negative integer
  1284.  */
  1285. static
  1286. int
  1287. hyperg_1F1_U(const double a, const double b, const double x, gsl_sf_result * result)
  1288. {
  1289.   const double bp = 2.0 - b;
  1290.   const double ap = a - b + 1.0;
  1291.  
  1292.   gsl_sf_result lg_ap, lg_bp;
  1293.   double sg_ap;
  1294.   int stat_lg0 = gsl_sf_lngamma_sgn_e(ap, &lg_ap, &sg_ap);
  1295.   int stat_lg1 = gsl_sf_lngamma_e(bp, &lg_bp);
  1296.   int stat_lg2 = GSL_ERROR_SELECT_2(stat_lg0, stat_lg1);
  1297.   double t1 = (bp-1.0) * log(x);
  1298.   double lnpre_val = lg_ap.val - lg_bp.val + t1;
  1299.   double lnpre_err = lg_ap.err + lg_bp.err + 2.0 * GSL_DBL_EPSILON * fabs(t1);
  1300.  
  1301.   gsl_sf_result lg_2mbp, lg_1papmbp;
  1302.   double sg_2mbp, sg_1papmbp;
  1303.   int stat_lg3 = gsl_sf_lngamma_sgn_e(2.0-bp,    &lg_2mbp,    &sg_2mbp);
  1304.   int stat_lg4 = gsl_sf_lngamma_sgn_e(1.0+ap-bp, &lg_1papmbp, &sg_1papmbp);
  1305.   int stat_lg5 = GSL_ERROR_SELECT_2(stat_lg3, stat_lg4);
  1306.   double lnc1_val = lg_2mbp.val - lg_1papmbp.val;
  1307.   double lnc1_err = lg_2mbp.err + lg_1papmbp.err
  1308.                     + GSL_DBL_EPSILON * (fabs(lg_2mbp.val) + fabs(lg_1papmbp.val));
  1309.  
  1310.   gsl_sf_result M;
  1311.   gsl_sf_result_e10 U;
  1312.   int stat_F = gsl_sf_hyperg_1F1_e(ap, bp, x, &M);
  1313.   int stat_U = gsl_sf_hyperg_U_e10_e(ap, bp, x, &U);
  1314.   int stat_FU = GSL_ERROR_SELECT_2(stat_F, stat_U);
  1315.  
  1316.   gsl_sf_result_e10 term_M;
  1317.   int stat_e0 = gsl_sf_exp_mult_err_e10_e(lnc1_val, lnc1_err,
  1318.                                              sg_2mbp*sg_1papmbp*M.val, M.err,
  1319.                                              &term_M);
  1320.  
  1321.   const double ombp = 1.0 - bp;
  1322.   const double Uee_val = U.e10*M_LN10;
  1323.   const double Uee_err = 2.0 * GSL_DBL_EPSILON * fabs(Uee_val);
  1324.   const double Mee_val = term_M.e10*M_LN10;
  1325.   const double Mee_err = 2.0 * GSL_DBL_EPSILON * fabs(Mee_val);
  1326.   int stat_e1;
  1327.  
  1328.   /* Do a little dance with the exponential prefactors
  1329.    * to avoid overflows in intermediate results.
  1330.    */
  1331.   if(Uee_val > Mee_val) {
  1332.     const double factorM_val = exp(Mee_val-Uee_val);
  1333.     const double factorM_err = 2.0 * GSL_DBL_EPSILON * (fabs(Mee_val-Uee_val)+1.0) * factorM_val;
  1334.     const double inner_val = term_M.val*factorM_val - ombp*U.val;
  1335.     const double inner_err =
  1336.         term_M.err*factorM_val + fabs(ombp) * U.err
  1337.       + fabs(term_M.val) * factorM_err
  1338.       + GSL_DBL_EPSILON * (fabs(term_M.val*factorM_val) + fabs(ombp*U.val));
  1339.     stat_e1 = gsl_sf_exp_mult_err_e(lnpre_val+Uee_val, lnpre_err+Uee_err,
  1340.                                        sg_ap*inner_val, inner_err,
  1341.                                        result);
  1342.   }
  1343.   else {
  1344.     const double factorU_val = exp(Uee_val - Mee_val);
  1345.     const double factorU_err = 2.0 * GSL_DBL_EPSILON * (fabs(Mee_val-Uee_val)+1.0) * factorU_val;
  1346.     const double inner_val = term_M.val - ombp*factorU_val*U.val;
  1347.     const double inner_err =
  1348.         term_M.err + fabs(ombp*factorU_val*U.err)
  1349.       + fabs(ombp*factorU_err*U.val)
  1350.       + GSL_DBL_EPSILON * (fabs(term_M.val) + fabs(ombp*factorU_val*U.val));
  1351.     stat_e1 = gsl_sf_exp_mult_err_e(lnpre_val+Mee_val, lnpre_err+Mee_err,
  1352.                                        sg_ap*inner_val, inner_err,
  1353.                                        result);
  1354.   }
  1355.  
  1356.   return GSL_ERROR_SELECT_5(stat_e1, stat_e0, stat_FU, stat_lg5, stat_lg2);
  1357. }
  1358.  
  1359.  
  1360. /* Handle case of generic positive a, b.
  1361.  * Assumes b-a is not a negative integer.
  1362.  */
  1363. static
  1364. int
  1365. hyperg_1F1_ab_pos(const double a, const double b,
  1366.                   const double x,
  1367.                   gsl_sf_result * result)
  1368. {
  1369.   const double ax = fabs(x);
  1370.  
  1371.   if(   ( b < 10.0 && a < 10.0 && ax < 5.0 )
  1372.      || ( b > a*ax )
  1373.      || ( b > a && ax < 5.0 )
  1374.     ) {
  1375.     return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  1376.   }
  1377.   else if(   x < -100.0
  1378.           && GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.7*fabs(x)
  1379.     ) {
  1380.     /* Large negative x asymptotic.
  1381.      */
  1382.     return hyperg_1F1_asymp_negx(a, b, x, result);
  1383.   }
  1384.   else if(   x > 100.0
  1385.           && GSL_MAX_DBL(fabs(b-a),1.0)*GSL_MAX_DBL(fabs(1.0-a),1.0) < 0.7*fabs(x)
  1386.     ) {
  1387.     /* Large positive x asymptotic.
  1388.      */
  1389.     return hyperg_1F1_asymp_posx(a, b, x, result);
  1390.   }
  1391.   else if(fabs(b-a) <= 1.0) {
  1392.     /* Directly handle b near a.
  1393.      */
  1394.     return hyperg_1F1_beps_bgt0(a-b, b, x, result);  /* a = b + eps */
  1395.   }
  1396.  
  1397.   else if(b > a && b >= 2*a + x) {
  1398.     /* Use the Gautschi CF series, then
  1399.      * recurse backward to a near 0 for normalization.
  1400.      * This will work for either sign of x.
  1401.      */ 
  1402.     double rap;
  1403.     int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
  1404.     double ra = 1.0 + x/a * rap;
  1405.  
  1406.     double Ma    = GSL_SQRT_DBL_MIN;
  1407.     double Map1 = ra * Ma;
  1408.     double Mnp1 = Map1;
  1409.     double Mn    = Ma;
  1410.     double Mnm1;
  1411.     gsl_sf_result Mn_true;
  1412.     int stat_Mt;
  1413.     double n;
  1414.     for(n=a; n>0.5; n -= 1.0) {
  1415.       Mnm1 = (n * Mnp1 - (2.0*n-b+x) * Mn) / (b-n);
  1416.       Mnp1 = Mn;
  1417.       Mn   = Mnm1;
  1418.     }
  1419.  
  1420.     stat_Mt = hyperg_1F1_small_a_bgt0(n, b, x, &Mn_true);
  1421.  
  1422.     result->val  = (Ma/Mn) * Mn_true.val;
  1423.     result->err  = fabs(Ma/Mn) * Mn_true.err;
  1424.     result->err += 2.0 * GSL_DBL_EPSILON * (fabs(a)+1.0) * fabs(result->val);
  1425.     return GSL_ERROR_SELECT_2(stat_Mt, stat_CF1);
  1426.   }
  1427.   else if(b > a && b < 2*a + x && b > x) {
  1428.     /* Use the Gautschi series representation of
  1429.      * the continued fraction. Then recurse forward
  1430.      * to near the a=b line for normalization. This will
  1431.      * work for either sign of x, although we do need
  1432.      * to check for b > x, which is relevant when x is positive.
  1433.      */
  1434.     gsl_sf_result Mn_true;
  1435.     int stat_Mt;
  1436.     double rap;
  1437.     int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
  1438.     double ra = 1.0 + x/a * rap;
  1439.     double Ma    = GSL_SQRT_DBL_MIN;
  1440.     double Mnm1 = Ma;
  1441.     double Mn    = ra * Mnm1;
  1442.     double Mnp1;
  1443.     double n;
  1444.     for(n=a+1.0; n<b-0.5; n += 1.0) {
  1445.       Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
  1446.       Mnm1 = Mn;
  1447.       Mn   = Mnp1;
  1448.     }
  1449.     stat_Mt = hyperg_1F1_beps_bgt0(n-b, b, x, &Mn_true);
  1450.     result->val  = Ma/Mn * Mn_true.val;
  1451.     result->err  = fabs(Ma/Mn) * Mn_true.err;
  1452.     result->err += 2.0 * GSL_DBL_EPSILON * (fabs(b-a)+1.0) * fabs(result->val);
  1453.     return GSL_ERROR_SELECT_2(stat_Mt, stat_CF1);
  1454.   }
  1455.   else if(x >= 0.0) {
  1456.  
  1457.     if(b < a) {
  1458.       /* Forward recursion on a from a=b+eps-1,b+eps.
  1459.        */
  1460.       double N   = floor(a-b);
  1461.       double eps = a - b - N;
  1462.       gsl_sf_result r_M0;
  1463.       gsl_sf_result r_M1;
  1464.       int stat_0 = hyperg_1F1_beps_bgt0(eps-1.0, b, x, &r_M0);
  1465.       int stat_1 = hyperg_1F1_beps_bgt0(eps,     b, x, &r_M1);
  1466.       double M0 = r_M0.val;
  1467.       double M1 = r_M1.val;
  1468.  
  1469.       double Mam1 = M0;
  1470.       double Ma   = M1;
  1471.       double Map1;
  1472.       double ap;
  1473.       double start_pair = fabs(M0) + fabs(M1);
  1474.       double minim_pair = GSL_DBL_MAX;
  1475.       double pair_ratio;
  1476.       double rat_0 = fabs(r_M0.err/r_M0.val);
  1477.       double rat_1 = fabs(r_M1.err/r_M1.val);
  1478.       for(ap=b+eps; ap<a-0.1; ap += 1.0) {
  1479.         Map1 = ((b-ap)*Mam1 + (2.0*ap-b+x)*Ma)/ap;
  1480.         Mam1 = Ma;
  1481.         Ma   = Map1;
  1482.         minim_pair = GSL_MIN_DBL(fabs(Mam1) + fabs(Ma), minim_pair);
  1483.       }
  1484.       pair_ratio = start_pair/minim_pair;
  1485.       result->val  = Ma;
  1486.       result->err  = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Ma);
  1487.       result->err += 2.0 * (rat_0 + rat_1) * pair_ratio*pair_ratio * fabs(Ma);
  1488.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(Ma);
  1489.       return GSL_ERROR_SELECT_2(stat_0, stat_1);
  1490.     }
  1491.     else {
  1492.       /* b > a
  1493.        * b < 2a + x 
  1494.        * b <= x
  1495.        *
  1496.        * Recurse forward on a from a=eps,eps+1.
  1497.        */
  1498.       double eps = a - floor(a);
  1499.       gsl_sf_result r_Mnm1;
  1500.       gsl_sf_result r_Mn;
  1501.       int stat_0 = hyperg_1F1_small_a_bgt0(eps,     b, x, &r_Mnm1);
  1502.       int stat_1 = hyperg_1F1_small_a_bgt0(eps+1.0, b, x, &r_Mn);
  1503.       double Mnm1 = r_Mnm1.val;
  1504.       double Mn   = r_Mn.val;
  1505.       double Mnp1;
  1506.  
  1507.       double n;
  1508.       double start_pair = fabs(Mn) + fabs(Mnm1);
  1509.       double minim_pair = GSL_DBL_MAX;
  1510.       double pair_ratio;
  1511.       double rat_0 = fabs(r_Mnm1.err/r_Mnm1.val);
  1512.       double rat_1 = fabs(r_Mn.err/r_Mn.val);
  1513.       for(n=eps+1.0; n<a-0.1; n++) {
  1514.         Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
  1515.         Mnm1 = Mn;
  1516.         Mn   = Mnp1;
  1517.         minim_pair = GSL_MIN_DBL(fabs(Mn) + fabs(Mnm1), minim_pair);
  1518.       }
  1519.       pair_ratio = start_pair/minim_pair;
  1520.       result->val  = Mn;
  1521.       result->err  = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(a)+1.0) * fabs(Mn);
  1522.       result->err += 2.0 * (rat_0 + rat_1) * pair_ratio*pair_ratio * fabs(Mn);
  1523.       result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn);
  1524.       return GSL_ERROR_SELECT_2(stat_0, stat_1);
  1525.     }
  1526.   }
  1527.   else {
  1528.     /* x < 0
  1529.      * b < a
  1530.      */
  1531.  
  1532.     if(a <= 0.5*(b-x) || a >= -x) {
  1533.       /* Recurse down in b, from near the a=b line, b=a+eps,a+eps-1.
  1534.        */
  1535.       double N   = floor(a - b);
  1536.       double eps = 1.0 + N - a + b;
  1537.       gsl_sf_result r_Manp1;
  1538.       gsl_sf_result r_Man;
  1539.       int stat_0 = hyperg_1F1_beps_bgt0(-eps,    a+eps,     x, &r_Manp1);
  1540.       int stat_1 = hyperg_1F1_beps_bgt0(1.0-eps, a+eps-1.0, x, &r_Man);
  1541.       double Manp1 = r_Manp1.val;
  1542.       double Man   = r_Man.val;
  1543.       double Manm1;
  1544.  
  1545.       double n;
  1546.       double start_pair = fabs(Manp1) + fabs(Man);
  1547.       double minim_pair = GSL_DBL_MAX;
  1548.       double pair_ratio;
  1549.       double rat_0 = fabs(r_Manp1.err/r_Manp1.val);
  1550.       double rat_1 = fabs(r_Man.err/r_Man.val);
  1551.       for(n=a+eps-1.0; n>b+0.1; n -= 1.0) {
  1552.         Manm1 = (-n*(1-n-x)*Man - x*(n-a)*Manp1)/(n*(n-1.0));
  1553.         Manp1 = Man;
  1554.         Man = Manm1;
  1555.         minim_pair = GSL_MIN_DBL(fabs(Manp1) + fabs(Man), minim_pair);
  1556.       }
  1557.  
  1558.       /* FIXME: this is a nasty little hack; there is some
  1559.          (transient?) instability in this recurrence for some
  1560.      values. I can tell when it happens, which is when
  1561.      this pair_ratio is large. But I do not know how to
  1562.      measure the error in terms of it. I guessed quadratic
  1563.      below, but it is probably worse than that.
  1564.      */
  1565.       pair_ratio = start_pair/minim_pair;
  1566.       result->val  = Man;
  1567.       result->err  = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Man);
  1568.       result->err *= pair_ratio*pair_ratio + 1.0;
  1569.       return GSL_ERROR_SELECT_2(stat_0, stat_1);
  1570.     }
  1571.     else {
  1572.       /* Pick a0 such that b ~= 2a0 + x, then
  1573.        * recurse down in b from a0,a0 to determine
  1574.        * the values near the line b=2a+x. Then recurse
  1575.        * forward on a from a0.
  1576.        */
  1577.       double epsa = a - floor(a);
  1578.       double a0   = floor(0.5*(b-x)) + epsa;
  1579.       double N    = floor(a0 - b);
  1580.       double epsb = 1.0 + N - a0 + b;
  1581.       double Ma0b;
  1582.       double Ma0bp1;
  1583.       double Ma0p1b;
  1584.       int stat_a0;
  1585.       double Mnm1;
  1586.       double Mn;
  1587.       double Mnp1;
  1588.       double n;
  1589.       double err_rat;
  1590.       {
  1591.         gsl_sf_result r_Ma0np1;
  1592.         gsl_sf_result r_Ma0n;
  1593.         int stat_0 = hyperg_1F1_beps_bgt0(-epsb,    a0+epsb,     x, &r_Ma0np1);
  1594.         int stat_1 = hyperg_1F1_beps_bgt0(1.0-epsb, a0+epsb-1.0, x, &r_Ma0n);
  1595.         double Ma0np1 = r_Ma0np1.val;
  1596.         double Ma0n   = r_Ma0n.val;
  1597.         double Ma0nm1;
  1598.  
  1599.     err_rat = fabs(r_Ma0np1.err/r_Ma0np1.val) + fabs(r_Ma0n.err/r_Ma0n.val);
  1600.  
  1601.         for(n=a0+epsb-1.0; n>b+0.1; n -= 1.0) {
  1602.           Ma0nm1 = (-n*(1-n-x)*Ma0n - x*(n-a0)*Ma0np1)/(n*(n-1.0));
  1603.           Ma0np1 = Ma0n;
  1604.           Ma0n = Ma0nm1;
  1605.         }
  1606.     Ma0bp1 = Ma0np1;
  1607.         Ma0b   = Ma0n;
  1608.     Ma0p1b = (b*(a0+x)*Ma0b+x*(a0-b)*Ma0bp1)/(a0*b); /* right-down hook */
  1609.     stat_a0 = GSL_ERROR_SELECT_2(stat_0, stat_1);
  1610.       }
  1611.  
  1612.           
  1613.       /* Initialise the recurrence correctly BJG */
  1614.  
  1615.       if (a0 >= a - 0.1)
  1616.         { 
  1617.           Mn = Ma0b;
  1618.         }
  1619.       else if (a0 + 1>= a - 0.1)
  1620.         {
  1621.           Mn = Ma0p1b;
  1622.         }
  1623.       else
  1624.         {
  1625.           Mnm1 = Ma0b;
  1626.           Mn   = Ma0p1b;
  1627.  
  1628.           for(n=a0+1.0; n<a-0.1; n += 1.0) {
  1629.             Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
  1630.             Mnm1 = Mn;
  1631.             Mn   = Mnp1;
  1632.           }
  1633.         }
  1634.  
  1635.       result->val = Mn;
  1636.       result->err = (err_rat + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Mn);
  1637.       return stat_a0;
  1638.     }
  1639.   }
  1640. }
  1641.  
  1642.  
  1643. /* Assumes b != integer
  1644.  * Assumes a != integer when x > 0
  1645.  * Assumes b-a != neg integer when x < 0
  1646.  */
  1647. static
  1648. int
  1649. hyperg_1F1_ab_neg(const double a, const double b, const double x,
  1650.                   gsl_sf_result * result)
  1651. {
  1652.   const double bma = b - a;
  1653.   const double abs_x = fabs(x);
  1654.   const double abs_a = fabs(a);
  1655.   const double abs_b = fabs(b);
  1656.   const double size_a = GSL_MAX(abs_a, 1.0);
  1657.   const double size_b = GSL_MAX(abs_b, 1.0);
  1658.   const int bma_integer = ( bma - floor(bma+0.5) < _1F1_INT_THRESHOLD );
  1659.  
  1660.   if(   (abs_a < 10.0 && abs_b < 10.0 && abs_x < 5.0)
  1661.      || (b > 0.8*GSL_MAX(fabs(a),1.0)*fabs(x))
  1662.     ) {
  1663.     return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  1664.   }
  1665.   else if(   x > 0.0
  1666.           && size_b > size_a
  1667.           && size_a*log(M_E*x/size_b) < GSL_LOG_DBL_EPSILON+7.0
  1668.     ) {
  1669.     /* Series terms are positive definite up until
  1670.      * there is a sign change. But by then the
  1671.      * terms are small due to the last condition.
  1672.      */
  1673.     return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  1674.   }
  1675.   else if(   (abs_x < 5.0 && fabs(bma) < 10.0 && abs_b < 10.0)
  1676.           || (b > 0.8*GSL_MAX_DBL(fabs(bma),1.0)*abs_x)
  1677.     ) {
  1678.     /* Use Kummer transformation to render series safe.
  1679.      */
  1680.     gsl_sf_result Kummer_1F1;
  1681.     int stat_K = gsl_sf_hyperg_1F1_series_e(bma, b, -x, &Kummer_1F1);
  1682.     int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
  1683.                                       Kummer_1F1.val, Kummer_1F1.err,
  1684.                                       result);
  1685.     return GSL_ERROR_SELECT_2(stat_e, stat_K);
  1686.   }
  1687.   else if(   x < -30.0
  1688.           && GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.99*fabs(x)
  1689.     ) {
  1690.     /* Large negative x asymptotic.
  1691.      * Note that we do not check if b-a is a negative integer.
  1692.      */
  1693.     return hyperg_1F1_asymp_negx(a, b, x, result);
  1694.   }
  1695.   else if(   x > 100.0
  1696.           && GSL_MAX_DBL(fabs(bma),1.0)*GSL_MAX_DBL(fabs(1.0-a),1.0) < 0.99*fabs(x)
  1697.     ) {
  1698.     /* Large positive x asymptotic.
  1699.      * Note that we do not check if a is a negative integer.
  1700.      */
  1701.     return hyperg_1F1_asymp_posx(a, b, x, result);
  1702.   }
  1703.   else if(x > 0.0 && !(bma_integer && bma > 0.0)) {
  1704.     return hyperg_1F1_U(a, b, x, result);
  1705.   }
  1706.   else {
  1707.     /* FIXME:  if all else fails, try the series... BJG */
  1708.     if (x < 0.0) {
  1709.       /* Apply Kummer Transformation */
  1710.       int status = gsl_sf_hyperg_1F1_series_e(b-a, b, -x, result);
  1711.       double K_factor = exp(x);
  1712.       result->val *= K_factor;
  1713.       result->err *= K_factor;
  1714.       return status;
  1715.     } else {
  1716.       int status = gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  1717.       return status;
  1718.     }
  1719.  
  1720.     /* Sadness... */
  1721.     /* result->val = 0.0; */
  1722.     /* result->err = 0.0; */
  1723.     /* GSL_ERROR ("error", GSL_EUNIMPL); */
  1724.   }
  1725. }
  1726.  
  1727.  
  1728. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  1729.  
  1730. int
  1731. gsl_sf_hyperg_1F1_int_e(const int a, const int b, const double x, gsl_sf_result * result)
  1732. {
  1733.   /* CHECK_POINTER(result) */
  1734.  
  1735.   if(x == 0.0) {
  1736.     result->val = 1.0;
  1737.     result->err = 0.0;
  1738.     return GSL_SUCCESS;
  1739.   }
  1740.   else if(a == b) {
  1741.     return gsl_sf_exp_e(x, result);
  1742.   }
  1743.   else if(b == 0) {
  1744.     DOMAIN_ERROR(result);
  1745.   }
  1746.   else if(a == 0) {
  1747.     result->val = 1.0;
  1748.     result->err = 0.0;
  1749.     return GSL_SUCCESS;
  1750.   }
  1751.   else if(b < 0 && (a < b || a > 0)) {
  1752.     /* Standard domain error due to singularity. */
  1753.     DOMAIN_ERROR(result);
  1754.   }
  1755.   else if(x > 100.0  && GSL_MAX_DBL(1.0,fabs(b-a))*GSL_MAX_DBL(1.0,fabs(1-a)) < 0.5 * x) {
  1756.     /* x -> +Inf asymptotic */
  1757.     return hyperg_1F1_asymp_posx(a, b, x, result);
  1758.   }
  1759.   else if(x < -100.0 && GSL_MAX_DBL(1.0,fabs(a))*GSL_MAX_DBL(1.0,fabs(1+a-b)) < 0.5 * fabs(x)) {
  1760.     /* x -> -Inf asymptotic */
  1761.     return hyperg_1F1_asymp_negx(a, b, x, result);
  1762.   }
  1763.   else if(a < 0 && b < 0) {
  1764.     return hyperg_1F1_ab_negint(a, b, x, result);
  1765.   }
  1766.   else if(a < 0 && b > 0) {
  1767.     /* Use Kummer to reduce it to the positive integer case.
  1768.      * Note that b > a, strictly, since we already trapped b = a.
  1769.      */
  1770.     gsl_sf_result Kummer_1F1;
  1771.     int stat_K = hyperg_1F1_ab_posint(b-a, b, -x, &Kummer_1F1);
  1772.     int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
  1773.                                       Kummer_1F1.val, Kummer_1F1.err,
  1774.                                       result); 
  1775.     return GSL_ERROR_SELECT_2(stat_e, stat_K);
  1776.   }
  1777.   else {
  1778.     /* a > 0 and b > 0 */
  1779.     return hyperg_1F1_ab_posint(a, b, x, result);
  1780.   }
  1781. }
  1782.  
  1783.  
  1784. int
  1785. gsl_sf_hyperg_1F1_e(const double a, const double b, const double x,
  1786.                        gsl_sf_result * result
  1787.                        )
  1788. {
  1789.   const double bma = b - a;
  1790.   const double rinta = floor(a + 0.5);
  1791.   const double rintb = floor(b + 0.5);
  1792.   const double rintbma = floor(bma + 0.5);
  1793.   const int a_integer   = ( fabs(a-rinta) < _1F1_INT_THRESHOLD && rinta > INT_MIN && rinta < INT_MAX );
  1794.   const int b_integer   = ( fabs(b-rintb) < _1F1_INT_THRESHOLD && rintb > INT_MIN && rintb < INT_MAX );
  1795.   const int bma_integer = ( fabs(bma-rintbma) < _1F1_INT_THRESHOLD && rintbma > INT_MIN && rintbma < INT_MAX );
  1796.   const int b_neg_integer   = ( b < -0.1 && b_integer );
  1797.   const int a_neg_integer   = ( a < -0.1 && a_integer );
  1798.   const int bma_neg_integer = ( bma < -0.1 &&  bma_integer );
  1799.  
  1800.   /* CHECK_POINTER(result) */
  1801.  
  1802.   if(x == 0.0) {
  1803.     /* Testing for this before testing a and b
  1804.      * is somewhat arbitrary. The result is that
  1805.      * we have 1F1(a,0,0) = 1.
  1806.      */
  1807.     result->val = 1.0;
  1808.     result->err = 0.0;
  1809.     return GSL_SUCCESS;
  1810.   }
  1811.   else if(b == 0.0) {
  1812.     DOMAIN_ERROR(result);
  1813.   }
  1814.   else if(a == 0.0) {
  1815.     result->val = 1.0;
  1816.     result->err = 0.0;
  1817.     return GSL_SUCCESS;
  1818.   }
  1819.   else if(a == b) {
  1820.     /* case: a==b; exp(x)
  1821.      * It's good to test exact equality now.
  1822.      * We also test approximate equality later.
  1823.      */
  1824.     return gsl_sf_exp_e(x, result);
  1825.   }
  1826.   else if(fabs(b) < _1F1_INT_THRESHOLD) {
  1827.     /* Note that neither a nor b is zero, since
  1828.      * we eliminated that with the above tests.
  1829.      */
  1830.     if(fabs(a) < _1F1_INT_THRESHOLD) {
  1831.       /* a and b near zero: 1 + a/b (exp(x)-1)
  1832.        */
  1833.       gsl_sf_result exm1;
  1834.       int stat_e = gsl_sf_expm1_e(x, &exm1);
  1835.       double sa = ( a > 0.0 ? 1.0 : -1.0 );
  1836.       double sb = ( b > 0.0 ? 1.0 : -1.0 );
  1837.       double lnab = log(fabs(a/b)); /* safe */
  1838.       gsl_sf_result hx;
  1839.       int stat_hx = gsl_sf_exp_mult_err_e(lnab, GSL_DBL_EPSILON * fabs(lnab),
  1840.                              sa * sb * exm1.val, exm1.err,
  1841.                              &hx);
  1842.       result->val = (hx.val == GSL_DBL_MAX ? hx.val : 1.0 + hx.val);  /* FIXME: excessive paranoia ? what is DBL_MAX+1 ?*/
  1843.       result->err = hx.err;
  1844.       return GSL_ERROR_SELECT_2(stat_hx, stat_e);
  1845.     }
  1846.     else {
  1847.       /* b near zero and a not near zero
  1848.        */
  1849.       const double m_arg = 1.0/(0.5*b);
  1850.       gsl_sf_result F_renorm;
  1851.       int stat_F = hyperg_1F1_renorm_b0(a, x, &F_renorm);
  1852.       int stat_m = gsl_sf_multiply_err_e(m_arg, 2.0 * GSL_DBL_EPSILON * m_arg,
  1853.                                             0.5*F_renorm.val, 0.5*F_renorm.err,
  1854.                             result);
  1855.       return GSL_ERROR_SELECT_2(stat_m, stat_F);
  1856.     }
  1857.   }
  1858.   else if(a_integer && b_integer) {
  1859.     /* Check for reduction to the integer case.
  1860.      * Relies on the arbitrary "near an integer" test.
  1861.      */
  1862.     return gsl_sf_hyperg_1F1_int_e((int)rinta, (int)rintb, x, result);
  1863.   }
  1864.   else if(b_neg_integer && !(a_neg_integer && a > b)) {
  1865.     /* Standard domain error due to
  1866.      * uncancelled singularity.
  1867.      */
  1868.     DOMAIN_ERROR(result);
  1869.   }
  1870.   else if(a_neg_integer) {
  1871.     return hyperg_1F1_a_negint_lag((int)rinta, b, x, result);
  1872.   }
  1873.   else if(b > 0.0) {
  1874.     if(-1.0 <= a && a <= 1.0) {
  1875.       /* Handle small a explicitly.
  1876.        */
  1877.       return hyperg_1F1_small_a_bgt0(a, b, x, result);
  1878.     }
  1879.     else if(bma_neg_integer) {
  1880.       /* Catch this now, to avoid problems in the
  1881.        * generic evaluation code.
  1882.        */
  1883.       gsl_sf_result Kummer_1F1;
  1884.       int stat_K = hyperg_1F1_a_negint_lag((int)rintbma, b, -x, &Kummer_1F1);
  1885.       int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
  1886.                             Kummer_1F1.val, Kummer_1F1.err,
  1887.                             result);
  1888.       return GSL_ERROR_SELECT_2(stat_e, stat_K);
  1889.     }
  1890.     else if(a < 0.0) {
  1891.       /* Use Kummer to reduce it to the generic positive case.
  1892.        * Note that b > a, strictly, since we already trapped b = a.
  1893.        * Also b-(b-a)=a, and a is not a negative integer here,
  1894.        * so the generic evaluation is safe.
  1895.        */
  1896.       gsl_sf_result Kummer_1F1;
  1897.       int stat_K = hyperg_1F1_ab_pos(b-a, b, -x, &Kummer_1F1);
  1898.       int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
  1899.                             Kummer_1F1.val, Kummer_1F1.err,
  1900.                             result);
  1901.       return GSL_ERROR_SELECT_2(stat_e, stat_K);
  1902.     }
  1903.     else {
  1904.       /* a > 0.0 */
  1905.       return hyperg_1F1_ab_pos(a, b, x, result);
  1906.     }
  1907.   }
  1908.   else {
  1909.     /* b < 0.0 */
  1910.  
  1911.     if(bma_neg_integer && x < 0.0) {
  1912.       /* Handle this now to prevent problems
  1913.        * in the generic evaluation.
  1914.        */
  1915.       gsl_sf_result K;
  1916.       int stat_K;
  1917.       int stat_e;
  1918.       if(a < 0.0) {
  1919.         /* Kummer transformed version of safe polynomial.
  1920.      * The condition a < 0 is equivalent to b < b-a,
  1921.      * which is the condition required for the series
  1922.      * to be positive definite here.
  1923.          */
  1924.         stat_K = hyperg_1F1_a_negint_poly((int)rintbma, b, -x, &K);
  1925.       }
  1926.       else {
  1927.         /* Generic eval for negative integer a. */
  1928.     stat_K = hyperg_1F1_a_negint_lag((int)rintbma, b, -x, &K);
  1929.       }
  1930.       stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
  1931.                                         K.val, K.err,
  1932.                                         result);
  1933.       return GSL_ERROR_SELECT_2(stat_e, stat_K);
  1934.     }
  1935.     else if(a > 0.0) {
  1936.       /* Use Kummer to reduce it to the generic negative case.
  1937.        */
  1938.       gsl_sf_result K;
  1939.       int stat_K = hyperg_1F1_ab_neg(b-a, b, -x, &K);
  1940.       int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
  1941.                             K.val, K.err,
  1942.                             result);
  1943.       return GSL_ERROR_SELECT_2(stat_e, stat_K);
  1944.     }
  1945.     else {
  1946.       return hyperg_1F1_ab_neg(a, b, x, result);
  1947.     }
  1948.   }
  1949. }
  1950.  
  1951.  
  1952.   
  1953. #if 0  
  1954.     /* Luke in the canonical case.
  1955.    */
  1956.   if(x < 0.0 && !a_neg_integer && !bma_neg_integer) {
  1957.     double prec;
  1958.     return hyperg_1F1_luke(a, b, x, result, &prec);
  1959.   }
  1960.  
  1961.  
  1962.   /* Luke with Kummer transformation.
  1963.    */
  1964.   if(x > 0.0 && !a_neg_integer && !bma_neg_integer) {
  1965.     double prec;
  1966.     double Kummer_1F1;
  1967.     double ex;
  1968.     int stat_F = hyperg_1F1_luke(b-a, b, -x, &Kummer_1F1, &prec);
  1969.     int stat_e = gsl_sf_exp_e(x, &ex);
  1970.     if(stat_F == GSL_SUCCESS && stat_e == GSL_SUCCESS) {
  1971.       double lnr = log(fabs(Kummer_1F1)) + x;
  1972.       if(lnr < GSL_LOG_DBL_MAX) {
  1973.         *result = ex * Kummer_1F1;
  1974.     return GSL_SUCCESS;
  1975.       }
  1976.       else {
  1977.         *result = GSL_POSINF; 
  1978.     GSL_ERROR ("overflow", GSL_EOVRFLW);
  1979.       }
  1980.     }
  1981.     else if(stat_F != GSL_SUCCESS) {
  1982.       *result = 0.0;
  1983.       return stat_F;
  1984.     }
  1985.     else {
  1986.       *result = 0.0;
  1987.       return stat_e;
  1988.     }
  1989.   }
  1990. #endif
  1991.  
  1992.  
  1993.  
  1994. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  1995.  
  1996. #include "eval.h"
  1997.  
  1998. double gsl_sf_hyperg_1F1_int(const int m, const int n, double x)
  1999. {
  2000.   EVAL_RESULT(gsl_sf_hyperg_1F1_int_e(m, n, x, &result));
  2001. }
  2002.  
  2003. double gsl_sf_hyperg_1F1(double a, double b, double x)
  2004. {
  2005.   EVAL_RESULT(gsl_sf_hyperg_1F1_e(a, b, x, &result));
  2006. }
  2007.